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A numerical scheme to study the mixed states in a mesoscopic type-II superconducting cylinder is 
described. Steady-state configurations and transient behavior of the magnetic vortices for various 
values of the applied magnetic field H are presented. Transitions between different multi-vortex 
states as H is changed is demonstrated by contour plots and jumps in the B vs H plot. Evolving 
a uniformly-superconducting initial state using the simplest set of relaxation equations shows that 
the system passes through nearly metastable intermediate configurations, while seeking the final 
minimum-energy, steady state consistent with the square symmetry of the sample. An efficient 
scheme to determine the equilibrium vortex configuration in a mesoscopic system at any given 
applied field, not limited to the symmetry of the system, is devised and demonstrated. 
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I. INTRODUCTION 

The electromagnetic characteristics of type-II super- 
conductors have been of great interest, particularly 
since the discovery of high-temperature superconduc- 
tors. Furthermore, a recent surge in interest about 
nano-technology, has drawn the attention of researchers 
from various disciplines to a detailed understanding of 
the characteristics of mesoscopic superconducting sam- 
ples. Several phenomenological theories have been de- 
veloped during decades of superconductor research, one 
popular choice is the Ginzburg-Landau theory.^'^ Its 
time-dependent extension is known as Time-Dependent 
Ginzburg-Landau (TDGL) theory.^ 

In the Ginzburg-Landau theory, the electromagnetic 
state of a superconductor can be determined by solving 
a system of partial differential equations. It was discov- 
ered by Abrikosov that if the k parameter that appears 
in the equations, now known as the Ginzburg-Landau 
parameter, is larger than l/-\/2, then when a bulk su- 
perconductor is placed in a sufficiently large magnetic 
field, the magnetic field penetrates the superconductor 
in the form of singly-quantized vortices. Around each 
vortex flows a supercurrent,^ confining a single quantum 
of magnetic flux within it. Superconductors with this 
property are known as type-II. 

Here we present the result of a numerical study 
about the magnetization process inside a superconduct- 
ing cylinder with sub-micron lateral dimension in an ex- 
ternal magnetic field. We have restricted the work to a 
square cross section of a linear size equal to 4.65 times A 
(the magnetic penetration depth). A typical value for A 



is about 500 A, and then the cross-sectional area is about 
0.054 fim'^. 

Previous works on the magnetization of a mesoscopic 
superconductor without pinning centers have been re- 
ported by Peeters et al. and others^'^", who have 
presented extensive calculations on the superconducting 
state in mesoscopic, type-I, superconducting thin films. 
In most cases they found transitions between giant vortex 
states of different circulation quantum numbers L, with 
some multi-vortex states occasionally appearing as ther- 
modynamically stable states, but mostly as metastable 
states. These predictions appear to have already received 
some level of experimental confirmation, although some 
discrepancies still exist^-'^. (Ref.^*^ have mainly compared 
the energy of a "3-2" vortex-antivortex molecule state 
with that of a single off-centered vortex state at L = 1, 
as both evolve to the equilibrium state of a single vortex 
at the center.) Misko et al.^^ have studied both type- 
I and type-II mesoscopic trianglular cylinders, and have 
shown that a vortex-antivortex molecule appears only if 
the sample is type-I. They considered only one field value 
at which L = 2 is favored, and did not consider vortex 
configurational transitions as the field changes. 

Our aim is to simulate how vortices enter and settle 
in stable arrangements when a mesoscopic type-II super- 
conductor of a given symmetry is first cooled below the 
critical temperature, and then an external magnetic field 
is applied. This is often termed zero field cooling (ZFC). 
We find that only vortex numbers and configurations con- 
sistent with the sample symmetry can appear in this case. 
It is known that there exist global minimum-energy vor- 
tex configurations with reduced symmetry, correspond- 
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ingly to final equilibrium states at general values of the 
applied field. To find these equilibrium states we de- 
veloped, and demonstrate here, an efficient numerical 
scheme. 

Our approach is to solve a set of simplified (and dis- 
cretized) TDGL equations, in which the coupling to the 
electric field is neglected, and the superconducting order 
parameter and the magnetic field are assumed to relax 
with the same time scale. These assumptions arc not 
physical, but are acceptable here, since we are only in- 
terested in obtaining the final steady-state vortex config- 
urations, and the symmetry-related qualitative behavior 
of the transient configurations and their evolution. (For a 
more physical set of TDGL equations see Tinkham"*. For 
an example of numerical solution of such a set of TDGL 
equations, see^^.) However, since the equations we have 
solved do not contain any thermal fluctuation terms, and 
the sample we considered has a perfect square symmetry, 
we find that when starting with the Meissner state with 
no field penetration, then the final steady-state vortex 
configurations we obtain all have perfect square symme- 
try, with vortex numbers also limited to only multiples 
of four (the symmetry number). These configurations 
would correspond to physical situations under zero-field 
cooling, if the physical sample has indeed perfect symme- 
try, and the temperature is sufficiently low, so that ther- 
mal fluctuations are not able to overcome any energy bar- 
rier for vortex entry, expulsion, or rearrangement. How- 
ever, if the sample surface has slight imperfection, or if 
the temperature is not sufficiently low, then these con- 
figurations arc, in most cases, not in equilibrium at the 
given magnetic field strength. Even the vortex number 
may not be correct. However, if we insert into the equa- 
tions terms to simulate thermal fluctuation, as in the 
method of simulated annealing, ^^'^^ then the simulation 
program will take a much longer time to run, and may 
become impractical even with a supercomputer. So we 
have devised an efficient scheme to find the equilibrium 
vortex configurations: Wo solve the same set of relax- 
ation equations without any thermal fluctuation terms. 
But instead of starting the solution with the Meissner 
state as the initial state, we devise artificial initial states 
with a given number of vortices in random positions. 
We present analytic expressions for such initial states, 
in terms of a widely known approximate expression for 
a singly-quantized vortex in cylinder coordinates. Then 
for vortex numbers not too different from the equilibrium 
number, the flnal steady states obtained by solving our 
relaxation equations will, in most cases, have the num- 
ber of vortices close to those of the initial states. By 
comparing the total Gibbs energies of these steady states 
with different vortex numbers — sometimes we obtain 
more than one configuration for the same vortex num- 
ber, (when the vortex number exceeds four,) then their 
Gibbs energies are also compared — we can find the state 
with the lowest total Gibbs energy, which we identify as 
the equilibrium state with the equilibrium vortex num- 
ber. We give an explicit demonstration of this scheme^^. 



that might be very useful in view of the recent inter- 
est in nanoscience and nanotechnology. In some nano- 
applications one may need to find the equilibrium vor- 
tex configurations in a mesoscopic superconducting sam- 
ple, at different applied magnetic fields. We note that 
in a bulk sample vortices like; to form a triangular lat- 
tice. Thus, when the sample does not conform with this 
symmetry, and if the sample is sufficiently small so that 
the boundary effect on the equilibrium vortex configura- 
tions is important, then the system is frustrated, and the 
equilibrium vortex configurations can be quite intriguing, 
and difficult to foresee. The scheme devised here is then 
very convenient and useful for finding the answer. 



II. TIME-DEPENDENT GINZBURG-LANDAU 
MODEL 

In an external magnetic field H, the Gibbs free energy 
density 5 of a superconducting state is given by^ 
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where /„ is the free energy density in the normal state 
in the absence of the field, \1/ is the complex-valued order 
parameter, with the superscript * denoting complex con- 
jugation, A the magnetic vector potential, h = V x A 
the induced magnetic field and H the applied mag- 
netic field. The supercurrent density is expressed as 

J^=VxVxA = ^ (**V* - *V**) - A. 
In the above is the " effective charge" of a Cooper pair 
which is twice the charge of an electron, and rUg its "ef- 
fective mass" which can be selected arbitrarily, but the 
conventional choice is twice the mass of an electron. Also, 
c is the speed of light, and Ti = h/2-K where h is Planck's 
constant. 

Ginzburg-Landau theory postulates that the Gibbs 
free energy of a superconducting sample O, G (^', A) = 
J^gdQ is at a minimum in the superconducting state. 
The celebrated Ginzburg-Landau equations are obtained 
by minimizing this functional with respect to ^ and A 
using the variational principle. 

Since a constant term does not change the end result 
of the variational technique, an algebraic manipulation is 
made to subtract /„ and add H • H/Stt to the g above. 
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Nondimensionalizing variables as 
/ X xj' H h' ~ ^ 
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The characteristic scales are: 
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I^'ol = \J—a/ (3 which is the magnitude of that mini- 
mizes the free energy in the absence of a field; the ther- 

modynamic critical field strength ifc = ( 47r|a| |\l/o| 1 , 

which divides the normal state and superconducting state 
regions in Type-I superconductor phase diagram; the 

/ 2 \l/2 

London penetration depth A = {-^^^^^) ! tlie coher- 

/ .2 \l/2 

ence length ^ = ( 2ms \ oi \ ) ' Ginzburg-Landau 
parameter k, = A/^. 

We obtain the dimensionless gauge-invariant free en- 
ergy functional (omitting primes for convenience), 
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The simplified TDGL model we employ to find solu- 
tions of the static GL equations may be viewed as a 
gradient flow with the energy functional. That is, the 
variation of (^', A) w.r.t. time should be in the oppo- 
site direction of the gradient of the energy functional, 
^ = ^ = -lis witli time, t, in units of 

the only relaxation time of the equations. The natural 
boundary conditions are: (1) the continuity of the paral- 
lel component of the magnetic field across the boundary 
surface: (V x A) x n = H x n, (for two-dimensional 
problems only, see below) and (2) the vanishing gauge- 
invariant normal derivative of ^': — zA) VE" • n = 0,^° 
with n denoting the outward surface normal. 



III. DISCRETIZATION AND CALCULATION 
PROCEDURE 



For long cylindrical samples, we need only solve a 2- 
D problem. We take A = {A{x,y), B{x,y),0) and H = 
(0, 0, H) where F = (V x A)^ = || - |f . 

Defining the link variables as below, the gauge in- 
variance is preserved in discretizing the Gibbs free en- 
ergy and the consequent TDGL equations. W{x, y) = 
exp [in J^' A{<;,y)d(;) , V{x,y) ~ exp (in B{x,r])dT]) 
Noting that \da:{W*'i>)\ = \{d:, - ikA) and 
\dy (V**)| = \{dy- 1KB) *|, we have 

G(*,A) =^(^-|M/|2 + i|^,|4+|VxA-: 
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We discretize the free energy functional on a staggered 
grid in O shown in Fig. 1. 21.22 -pj^-g gj^gg second- 
order approximation in and hy 
ergy functional, where the hx and h 
crements in the x- and y-direction. The staggered grid 



to the continuous en- 
y are the spatial in- 
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FIG. 1. Staggered Grid 



also leads to a satisfactory way of discretizing the natural 
boundary conditions. For a rectangular grid, the first 
component of the vector potential is constant in time on 
one pair of the edges of the boundary, and the second 
component is constant in time on the other pair.^*^ 

In this paper, we assume that the cylindrical super- 
conductor has a square cross-section and is subject to an 
applied field along the cylindrical axis. The applied field 
is assumed to be constant in time. We further assume 
the order parameter ^E* varies in the cross-sectional plane 
of the cylindrical sample, and the vector potential A has 
only two nonzero components {A,B), which also lie in 
this plane. We also assume that the superconductor has 
no pinning sites. Then at steady-state conditions, the 
vortices settle in maximal distances due to the repulsion 
of each other. This requirement leads to triangular lat- 
tice of vortices in an infinitely large domain.^ 

In the staggered grid the lattice points for ^, A, and 
B are all different, with ^ evaluated at the node center 
(i,j), A evaluated at the east cell face (i-|-l/2, j), and B 
evaluated at the north cell face (i,j+l/2). According to 
Refs. 21, 22, this formulation keeps second order accuracy 
in the derivative evaluations as they appear in each of the 
discretized equations. 

Now the discrete TDGL is obtained by minimizing the 
discrete energy functional Gd with respect to the varia- 
tion in ^ and A as 
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BnE — + BnW . 

ft. ' 

— iV3(S„,*p,*A,) 
K 



Nx{^p) = (1- |*p|^) *p 



(7) 



(8) 



(^e, *p, = ($p0£; - Op^b) cos {A^Kh^) - 

($P$B + QpQe) sin (AeK/lcc) 

iVa *p, *Ar) = (^pOjv - Qp^n) cos (B^Khy) - 
($P$jv + 0p0jv) sin (BnKhy) 

where O and $ are the real and imaginary parts of ^, 

and boundary conditions of the computational domain D, 
(with T, B, L, and R denoting top, bottom, left, and 
right, respectively): 
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The finite difference equations arc solved by the Euler 
method with hx = hy = 0.15 and At = 0.05, and taking 
K = 4. In the numerical computations that follow, ev- 
erything is kept the same except for the strength of the 
applied magnetic field and/or the initial conditions. 



IV. STEADY STATES UNDER ZERO-FIELD 
COOLING IN A PERFECTLY SQUARE SAMPLE 
AT LOW TEMPERATURES 

We first solve the above set of equations assuming that 
the initial state is the perfect Meissner state with no field 
penetration. As explained in the introduction, this corre- 
sponds to applying a magnetic field after zero-field cool- 
ing. Fig. 2 shows the plots of l^fp [the left figure in 
(a) through (k)] and /i = V x A [the right figure in (a) 
through (k)] for the final steady states reached at a se- 
quence of increasing H values. In the left figures \^\^ , 
which is interpreted physically as the density of Cooper 



pairs, runs from to 1, the level 1 corresponding to the 
full superconducting state. Each isolated group of con- 
tours is called a "vortex" representing the supercurrent J 
circling around the vortex core, with ^ = Q at the vortex 
core. In the 3-D plots of Fig. 2 (b), it is clear that the 
vortices reach close to jvPl^ = at the core. 

Note that the number of vortices increases in multiples 
of 4. This is a consequence of the fact that the vortices 
are symmetrically created at the mid-points of the sample 
edges. Perfect symmetry in the sample geometry dictates 
that each side creates an equal number of vortices. The 
symmetry manifests itself in this geometry-dominated 
problem, and vortices arrange themselves in the square- 
symmetric configurations. The resultant steady states 
are mostly not true equilibrium states since the vortex 
arrangements do not reflect the intrinsic tendency of vor- 
tices to form a triangular lattice known to appear in bulk 
samples. The natural next step is to add a thermal fluctu- 
ation term to find the true equilibrium states which may 
or may not conform with this symmetry. This approach 
would then be like simulated annealing. ^^'^^ However, we 
believe it would not be practical to perfect this approach 
since it will likely be difficult to determine the appropri- 
ate rate of cooling and starting temperature. The nm- 
time of the computer program might also be expected to 
be much longer than we have found here. So we have de- 
vised a different approach which we believe is much more 
efficient for finding the equilibrium states. This is given 
in a later section. We shall see that even the cases with 
low number of vortices are not the true equilibrium. Also 
of interest is the fact that for some H the vortex config- 
uration has a much more difficult time to settle down, 
needing much longer run-times to get to a steady state. 
Geometry controls the settling time more than the energy 
in these cases. 

Our results arc summarized in Table 1, which lists the 
range of H for each possible number of vortices and the 
induced magnetic field -B = /fj hd^. The H's listed 
correspond to the threshold values for each range. They 
were found on a trial-and-error basis, and can be refined 
to any desired accuracy. Below H = 0.839, there is no 
vortex. Between H — 0.84 and 1.144, the vortex number 
n = 4, and so on. Of course, these thresholds change as 
we change the sample size. For example, if we double the 
length of each side, H = 1.0 gives n = 36. 
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H 


B 





0.839 


0.616350 


4 


0.84 


0.733198 


4 


1.144 


0.977913 


8 


1.145 


1.069432 


8 


1.429 


1.302241 


12 


1.43 


1.379098 


12 


1.732 


1.628862 


16 


1.733 


1.698629 


16 


2.058 


1.971973 



Table 1. The number of vortices n and the induced 
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(a) H =0.839 





(b) H =0.840 




(c) H =1.144 
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(d) H = 1.145 




(e) H = 1.429 




(f) H = 1 .430 



(g) H = 1 .454 
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(h) H = 1 .455 





(i) H = 1 .732 
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(j) H = 1 .733 




(k) H = 2.058 

FIG. 2. Plots of [left figures in (a) through (k)] and 
ft = V X A [right figures in (a) through (k)] for various H. 



magnetic field B for the applied magnetic field H. 

The B MS H plot is shown in Fig. 3, and reveals an 
abrupt increase in B between the regions of different 
number of vortices. For example, when H changes from 
1.144 to 1.145, B changes from 0.977913 to 1.069432. If 
the limit Ai? ^ is taken, we expect a sudden config- 
urational phase transition increasing the number of vor- 
tices, as is apparent in Figs. 2 (c) and (d). Such mini- 
first-order transitions are known to occur in a mesoscopic 
superconductor as H is changed, but the details are 
quite different, because different parameter (k) regimes 
and sample geometries (cylinder vs. film) are studied. 
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FIG. 3. plotted is B as a function of H . 

Phase transition is also evident in the n — 12 case, 
where the vortex configuration shows a sudden change in 
arrangement even with the same number of vortices for 
a slight increase of the applied field from H = 1.454 to 
H = 1.455. 

In the B vs H plot, B is much lower than H when the 
number of vortices is small, but as the vortices increase, 
the curve approaches the B = H curve asymptotically. 
So for the given sample, the number of vortices and B 
increase as the applied magnetic field H increases. This 
refiects the fact that the vortices are a form of magnetic 
penetration. 



V. TIME SEQUENCE SHOWING VORTEX 
ENTRY DYNAMICS 




FIG. 4. Time sequence of vortex-entry dynamics 
H = 1.145. 



for 



a square-symmetric configuration. We note that during 
the rotation process the vortex configuration loses some 
mirror symmetries of the sample but still preserves the 
90° rotation symmetry. So these transient states do not 
possess the full symmetry of the sample. We think that 
this is possible because our numerical method has very 
weekly broken the sample symmetry. That is, we think 
that the state just before the rotation is a metastable 
state only within the subspace of configurations preserv- 
ing the full symmetry of the sample. Thus, in the actual 
physical situation when the sample has perfect symmetry 
and the temperature is sufficiently low this rotation may 
take a very long time to take place. For samples with 
imperfect symmetry this relaxation time may be shorter. 
Since this is a symmetry-induced qualitative property of 
the vortex-entry dynamics in a mesoscopic superconduc- 
tor, we believe its general validity is independent of the 
fact that we have obtained it by solving a simplified set 
of TDGL equations which are not truly physical. 



Fig. 4 shows plots of Cooper pair density for H — 1.145 
in time. The number of vortices is 8. The perfect symme- 
try in the sample geometry dominates the transient pro- 
cess, but in the middle of the process the whole configu- 
ration makes a rotation to rearrange itself into a new con- 
figuration. (Note that time advances from 3000 to 17500 
between the 8th and 9th frames.) The final result is still 



VI. STEADY STATES WITH REDUCED 
SYMMETRY AND THE EQUILIBRIUM STATE 

The previous sections present solutions for a meso- 
scopic type-II superconducting cylinder with initially no 
vortex inside the system. The validity of such solutions 
also requires a perfectly square sample without any de- 
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feet at the boundary, and temperature sufficiently low, so 
that thermal fluctuations are too weak to help the sys- 
tem find lower-energy configurations of reduced symme- 
try. This is an ideal condition, producing only solutions 
consistent with the sample symmetry, and, even during 
the transient, the system is bound to this symmetry (ex- 
cept in rare cases when the transient solutions can keep 
only fourfold rotation symmetry but not mirror symme- 
tries See Fig. 4). In principle, one can reproduce this 
ideal system in a laboratory with special care. 

In real situations there most-likely exist some small 
defects or perturbations at the boundary. Then vortices 
can enter the system asymmetrically to produce steady- 
state configurations with reduced symmetry of lower to- 
tal Gibbs energy than any symmetric solution. A strong 
enough thermal fluctuation could also change the vortex 
number and rearrange the vortices to such a configura- 
tion. In previous work, to take into account these per- 
turbations, a random fluctuation term was added to the 
governing equation.^'* This term breaks the symmetry 
governing the equations, energizing the system to jump 
out of the local minima in energy, and over the energy 
barrier."^ But this increases the computing time greatly. 
(This method is basically what is called "simulated an- 
nealing" ^^'^^ in computer science.) 

As an alternative approach, we employ perturbed ini- 
tial conditions instead of the perfectly superconducting 
initial condition, and is similar to Peeters et al.^'"'' but 
with a new idea introduced to make the numerical scheme 
much more efficient: 

We have first used randomly perturbed initial condi- 
tions. They can indeed lead to final steady-state solu- 
tions with reduced symmetry and lower energies. But we 
find this way is very inefficient for finding the equilib- 
rium state at any given H . We have also tried to use a 
lower-symmetry configuration from such a calculation as 
the initial condition for a new H value, but find that the 
vortex number can often be trapped in an uncontrollable 
non-equilibrium value due to the existence of surface en- 
ergy barriers against vortex entry or exit. So this method 
of adopting an existing solution as the initial condition 
can not be reliably used to find the true equilibrium state 
in a given system and field. (Peeters et al. changes the 
field in small steps to avoid this diflaculty^^, but we be- 
lieve the procedure to be unnecessarily tedious.) 

To obtain the true equilibrium vortex configuration at 
any given magnetic field without employing some sim- 
ulated annealing method, we have devised a systematic 
approach to generate initial states with given numbers of 
vortices at random distributions. It is through an ana- 
lytic expression as follows: First, for one vortex at the 
origin in circular coordinates (r, 9) , we use the widely- 
known approximate expression^^"^^ 

re*" 

^{r,9) = ^===. (17) 

Converting it to Cartesian coordinates, we can move the 
center of the vortex to any arbitrary position {x',y') by 



simply replacing (x, y) by (x — x' ,y ~ y'). Denoting this 
expression as ^ x' ,v'{x,y), an n-vortex expression can be 
simply constructed as 

This expression obeys the important topological condi- 
tion that the phase of ^ must increase by 2tt when any 
one vortex center is circumnavigated. The magnetic field 
inside the sample does not obey any topological condi- 
tion, so it can be simply set equal to zero for the initial 
condition. The positions of the vortices can be gener- 
ated using random number generators, only if they are 
inside the sample. This idea may seem to be simple, but 
it appears to have not been employed before. We illus- 
trate below such initial conditions to obtain steady-state 
vortex configurations of any given numbers of vortices n. 
Comparing the total Gibbs energies of such solutions of 
different n, we can then determine the equilibrium vor- 
tex number and configuration. For illustrative purposes, 
we consider the case H — 0.840. In Fig. 5, the initial 
conditions for with 1 through 8 randomly-placed (ar- 
tificial) vortices [the left figures in (a) through (o)], and 
the steady-state vortex configurations they evolve to [the 
right figures in (a) through (o)], are shown. The corre- 
sponding Gibbs free energies of these steady states are 
plotted in Fig. 6 as a function of the vortex number n. 
The minimum-energy configuration at n = 5 is seen to 
display the square-symmetry of a five-vortex configura- 
tion with a vortex in the center. Although we have not 
yet applied this scheme to other field values, the method 
we have devised to find the equilibrium vortex config- 
urations for a given size and shape of the sample and 
different values of the external magnetic field should now 
be clear. 



VII. SUMMARY AND CONCLUSION 

A numerical scheme to study the mixed states in a 
mesoscopic type-II superconducting cylinder in a 

longitudinal external magnetic field H has been de- 
veloped. It is based on solving a set of simplified time- 
dependent Ginzburg-Landau equations. We have first 
applied this scheme to the case of field penetration into 
a zero-field cooled sample. Case studies for various val- 
ues of the external magnetic field are presented. Contour 
plots of the Cooper pair density, and the induced mag- 
netic field inside the sample, display the magnetic vortex 
solution first discovered by Abrikosov, but in a small sam- 
ple the vortex arrangement is not simply triangular. Gi- 
ant vortices and anti- vortices are not found in this study, 
unlike previous studies of type-I mesoscopic thin films. 
(But at sufficiently high magnetic field we still expect 
the system to favor a single giant vortex at the center 
as it goes into a surface superconducting state, but only 
if the sample is not too small.) Since we start the solu- 
tion with a uniformly superconducting initial condition, 
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OqO 



(n) 





(o) 

FIG. 5. The initial, random vortex configurations [left fig- 
ure in (a) through (o)], and the corresponding steady-state 
vortex configurations they evolve to [right figure in (a) 
through (o)]. 



-0.00725, 




and the sample has perfect square symmetry, both the 
number of vortices and their steady-state configurations 
are governed by the square sample geometry. Changes in 
the configuration and the number of vortices occur as H 
is varied through first-order configurational phase tran- 
sitions, similar to those found earlier, but different in de- 
tail. This phase transition characteristic is confirmed by 
the contour plots, and jumps in the values of the induced 
magnetic field B at certain discrete H values. A time se- 
quence shows that the system passes through intermedi- 
ate configurations, and remains in some of them for a long 
time, eventually settling down to the steady-state config- 
uration, which corresponds to the lowest- Gibbs-energy 
configuration consistent with the symmetry constraints 
to the vortex number and configuration. One could find 
the true equilibrium states, (which would appear in ac- 
tual samples, when there are symmetry-breaking sur- 
face defects as vortex-nucleation centers, or when ther- 
mal fluctuation is sufficiently strong to move the system 
out of metastable states, but not too strong to melt the 
vortex lattice,) by adding additional terms in the equa- 
tions to simulate thermal fluctuations, but here we have 
devised a different approach which we believe is more 
efficient. We introduce a way to generate analytic ini- 
tial states of prescribed numbers of vortices, but allow 
their positions to be random. They evolve to steady- 
state vortex arrangements of all possible vortex numbers 
near the equilibrium number, from which wc can compare 
total Gibbs energy to determine the equilibrium vortex 
number and configuration. In this way we can avoid the 
problem of surface and bulk energy barriers, which can 
trap the system in non-equilibrium vortex numbers and 
configurations — an undesirable situation which usually 
happens if one chooses the initial state randomly without 
controlling the vorticity quantum number L. 
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FIG. 6. The steady-state Gibbs free energy for various n. 
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